#!/bin/bash
set -e

while getopts  "p:" opts
do
        case  $opts  in
        
        
		p)
                        out_prefix=$OPTARG
                        
                ;;
		\?)

				echo $0 [-p out_prefix]
				;;
        esac
done
shift $(($OPTIND - 1))



if [ -z $out_prefix ]; then
	out_prefix=1
fi

#-----------------------------------------------
#-----------------------------------------------
. /mnt/ilustre/app/medical/tools/.var #---------
#-----------------------------------------------
#-----------------------------------------------

log=.log
if [ ! -e "$log" ]; then
	:> $log
fi

echo out prefix is: $out_prefix 2>>$log 1>&2


###--------------- argument may be changed ---------------###


# genome_name=b37.fa
echo using $genome_name as the reference genome 2>>$log 1>&2
# genome_assembly=b37
echo genome assembly is: $genome_assembly 2>>$log 1>&2
# dbsnp_version=138
echo dbsnp version is: $dbsnp_version 2>>$log 1>&2

# data_thread_num=8
# cpu_thread_num=4

# java_memory=16g
echo java memory: $java_memory 2>>$log 1>&2

# snpeff_db_version=GRCh37.75
echo snpeff database: $snpeff_db_version 2>>$log 1>&2

# header=\
# @HD\tVN:1.4\tGO:none\tSO:coordinate

# read_group=\
# @RG\\tID:${sample_name}\\tPL:ILLUMINA\\tSM:$sample_name
# For gatk, Each read group must contain the platform (PL) and sample (SM) tags.
# For the platform value, we currently support 454, LS454, Illumina, Solid, ABI_Solid, and CG (all case-insensitive).
# Each read in the file must be associated with exactly one read group.


vcf_path=${data_path}/vcf/gatk/

# there is no space between "=" and "\"
# ref_genome=\
# ${data_path}/ref/b37/$genome_name

echo ref genome is: $ref_genome 2>>$log 1>&2


reads_seq_1=$1
reads_seq_2=$2

echo reads are: $reads_seq_1 and $reads_seq_2 2>>$log 1>&2


###--------------- end argument may be changed ---------------###

###--------------- tools path ---------------###

# bwa=${tools_path}/bwa-0.7.12/bwa
# bwa=/share/apps/bwa-0.7.5a/bwa
# samtools=${tools_path}/samtools-1.2/samtools
# picard_path=${tools_path}/picard-tools-1.119/
# gatk_path=${tools_path}/GenomeAnalysisTK-1.6-9-g47df7bb/
# gatk3_path=${tools_path}
# gatk=${gatk3_path}/GenomeAnalysisTK.jar
# snpeff_path=${tools_path}/snpEff
# snpeff=${snpeff_path}/snpEff.jar
# snpsift=${snpeff_path}/SnpSift.jar
# fastuniq=${tools_path}/FastUniq/fastuniq
###--------------- end tools path ---------------###






echo 2>>$log 1>&2
echo 2>>$log 1>&2
echo gatk RealignerTargetCreator 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-T RealignerTargetCreator \
	-R $ref_genome \
	-I $out_prefix.sort.bam \
	-o $out_prefix.realn.intervals \
	-log $out_prefix.realn.log \
	-nt $data_thread_num \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	2>> $log
	
	
echo 2>>$log 1>&2
echo 2>>$log 1>&2
# indel realign
echo gatk IndelRealigner 2>>$log 1>&2
java -Xmx$java_memory -jar $gatk \
	-T IndelRealigner \
	-R $ref_genome \
	-targetIntervals $out_prefix.realn.intervals \
	-I $out_prefix.sort.bam \
	-o $out_prefix.realn.bam \
	-log $out_prefix.realn2.log \
	-known ${vcf_path}1000G_phase1.indels.${genome_assembly}.vcf \
	-known ${vcf_path}Mills_and_1000G_gold_standard.indels.${genome_assembly}.vcf \
	2>> $log